Motivation

Provide in one place some code to simulate data and fit spatial capture-recapture models, for both closed and open populations.

The reference book for spatial capture-capture (SCR) models is .

Slides of an introduction to SCR models by Andy Royle,

We will go the Bayesian way. Richard Chandler provides nice slides for Bayesian inference in closed population SCR models, and open population SCR models.

We will use Nimble to fit SCR models. We might implement a few tricks to improve and reach convergence faster.

Load Nimble.

library(nimble)

Closed population spatial capture-recapture models

Description

We use some code shared by Jose Jimenez on the SCR Google group and adapted by Perry de Valpine.

Simulations

The traps.

tr <- seq(15, 85, length = 10)
traps <- data.frame(X = rep(tr, each = length(tr)),
                    Y = rep(tr, times = length(tr))) # 100 coord. traps

Visualize.

viz_traps <- traps %>% 
  ggplot(aes(x = X, y = Y)) +
  geom_point(pch = 3) + 
  xlim(0, 100) +
  ylim(0, 100)
viz_traps

Generate population.

set.seed(10)
xlim <- c(0, 100)
ylim <- c(0, 100) # area 100 * 100 = 1e4
A <- (xlim[2] - xlim[1]) * (ylim[2] - ylim[1])/10000
A
## [1] 1
mu <- 50 # density
N <- rpois(1, mu*A) # generate population
N 
## [1] 50

Generate activity centers.

s <- data.frame(s.x = runif(N, xlim[1], xlim[2]), 
                s.y = runif(N, ylim[1], ylim[2]))

Visualize.

viz_traps_ac <- viz_traps + 
  geom_point(data = s, aes(x = s.x, y = s.y), pch = 16, color = "red")
viz_traps_ac

Generate detections.

sigma <- 5
lambda0 <- 0.4
J <- nrow(traps) # nb of traps
K <- 5 # nb capture occasions
yy <- array(NA, c(N, J, K))
for(j in 1:J) {
  dist <- sqrt((traps$X[j] - s$s.x)^2 + (traps$Y[j] - s$s.y)^2)
  lambda <- lambda0 * exp(-dist^2 / (2 * sigma^2))
  for(k in 1:K) {
    yy[,j,k] <- rpois(N, lambda)
  }
}
n <- apply(yy, c(2,3), sum)

Plot detections.

tot <- apply(n, 1, sum)
dat <- data.frame(traps, tot = tot)

viz_traps_ac +
  geom_point(data = dat, aes(x = X, y = Y, size = tot), alpha = 0.3) +
  scale_size(range = c(0, 20)) +
  labs(x = "",
       y = "",
       size = "# detections")

Model fitting

Define the model.

code <- nimbleCode({
  sigma ~ dunif(0, 10)
  lam0 ~ dunif(0, 5)
  psi ~ dbeta(1, 1)
  for(i in 1:M) {
    z[i] ~ dbern(psi)
    s[i,1] ~ dunif(xlim[1], xlim[2])
    s[i,2] ~ dunif(ylim[1], ylim[2])
    dist[i,1:J] <- (s[i,1] - X[1:J,1])^2 + (s[i,2] - X[1:J,2])^2
    lam[i,1:J] <- exp(-dist[i,1:J] / (2 * sigma^2)) * z[i]
  }
  for(j in 1:J){
    bigLambda[j] <- lam0 * sum(lam[1:M,j])
    for(k in 1:K) {
      n[j,k] ~ dpois(bigLambda[j])
    }
    }
  N <- sum(z[1:M])
})

Define constants, data and inits.

M <- 200
constants <- list(M = M, 
                  K = K, 
                  J = J)
n1 <- apply(n, 1, sum)
data <- list(n = n, 
             X = traps, 
             xlim = xlim, 
             ylim = ylim)
s <- cbind(runif(M, xlim[1], xlim[2]), 
           runif(M, ylim[1], ylim[2]))
z <- rep(1, M)
inits <- list(sigma = 0.5, 
              lam0 = 0.1, 
              s = s, 
              z = z,
              psi = 0.5)

Build R model (not compiled yet).

Rmodel <- nimbleModel(code = code, 
                      constants = constants, 
                      data = data, 
                      inits = inits)

Check whether the model is fully initialized. If you failed at providing initial values for some parameters (e.g. \(\psi\)), you’ll get NAs.

Rmodel$calculate()
## [1] -8089.009

Now compile the model in C++.

Cmodel <- compileNimble(Rmodel)

The R and C models are exactly the same versions of the model.

calculate(Cmodel)
## [1] -8089.009

You can simulate from prior.

Cmodel$simulate('lam0')
calculate(Cmodel)
## [1] -7669.788

Specify MCMC.

conf <- configureMCMC(Rmodel,
                      monitors = c("N", "lam0", "psi", "sigma"))
## ===== Monitors =====
## thin = 1: N, lam0, psi, sigma
## ===== Samplers =====
## RW sampler (402)
##   - sigma
##   - lam0
##   - s[]  (400 elements)
## conjugate sampler (1)
##   - psi
## binary sampler (200)
##   - z[]  (200 elements)

Build an executable MCMC.

Rmcmc <- buildMCMC(conf)

Compile in C++.

Cmcmc <- compileNimble(Rmcmc, project = Cmodel)

Run compiled model (do not run th uncompiled model with runMCMC(Rmcmc,100)).

samples <- runMCMC(Cmcmc,100)
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|

Explore.

dim(samples)
## [1] 100   4
colnames(samples)
## [1] "N"     "lam0"  "psi"   "sigma"
samplesSummary(samples)
##              Mean      Median    St.Dev.  95%CI_low   95%CI_upp
## N     126.4300000 123.5000000 19.7638866 95.9500000 185.3500000
## lam0    1.0657507   1.0453553  0.1245833  0.7794961   1.2568560
## psi     0.6426589   0.6284448  0.1080914  0.4638732   0.9672855
## sigma   2.2640864   2.2800982  0.1605662  1.7889305   2.4354455

Run compiled model.

samplesList <- runMCMC(Cmcmc,
                   niter = 10000,
                   nburnin = 5000,
                   nchains = 2)
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
## |-------------|-------------|-------------|-------------|
## |-------------------------------------------------------|
samples <- rbind(samplesList[[1]],
                 samplesList[[2]])
str(samples)
##  num [1:10000, 1:4] 31 29 27 23 22 22 19 24 28 29 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:4] "N" "lam0" "psi" "sigma"

Calculate ESS effective sample size (should be at least a thousand). Efficiency measures the tradeoff between ESS and computation time (the ratio). For example, an efficiency of 2 means that for every second you run the algorithm, it produces 2 independent samples.

library(coda)
apply(samples, 2, effectiveSize)
##         N      lam0       psi     sigma 
##  61.64045 143.41820  58.91601  44.60710

Produce trace an density plots.

library(basicMCMCplots)
chainsPlot(samplesList,
           var = c("N", "sigma", "lam0"))

Display summary stats. Compare to the values used to simulate data, in particuler \(N = 50\), \(\sigma = 5\) and \(\lambda_0 = 0.4\).

summary(samples)
##        N               lam0             psi              sigma      
##  Min.   :  9.00   Min.   :0.2070   Min.   :0.02896   Min.   :2.773  
##  1st Qu.: 38.00   1st Qu.:0.4230   1st Qu.:0.19116   1st Qu.:3.965  
##  Median : 55.00   Median :0.5031   Median :0.27948   Median :4.504  
##  Mean   : 61.08   Mean   :0.5151   Mean   :0.30760   Mean   :4.619  
##  3rd Qu.: 80.00   3rd Qu.:0.5903   3rd Qu.:0.40230   3rd Qu.:5.076  
##  Max.   :174.00   Max.   :1.3284   Max.   :0.87699   Max.   :9.349

Open population spatial capture-recapture models

Description

Simulations

Model fitting